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The problem of the evolution of a heavy quarkonium in a medium can be recast as that of 
a quantum dissipative system. Within the framework of the master-equation approach to open 
quantum systems, we consider the real-time dynamics of quarkonia. We find that in a plasma at 
fixed temperature, the populations of the various quarkonium states evolve together, while their 
momentum distribution satisfies a Fokker-Planck equation. 
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I. INTRODUCTION 

The original idea that heavy quarkonia might be sup- 
pressed in a deconfmed QCD medium, thereby probing 
the formation of a quark-gluon plasma in high-energy 
nuclear collisions [lj, has motivated many studies (see 
Refs. [2|-|4| for recent reviews). The properties of quarko- 
nia in a medium — be it deconfined or hadronic, as it was 
realized that the most fragile states might already be 
destroyed in a hot hadron gas — have been extensively in- 
vestigated, both in lattice QCD calculations and using 
effective field-theoretical approaches j5[. 

In particular, it was understood that a description of 
heavy quarkonia in a medium as non-relativistic quark- 
antiquark systems bound by an effective static poten- 
tial is possible, provided the potential has an imaginary 
part p-QjJ , which accounts for the finite lifetime of the 
states. This emphasizes the necessity to consider the 
real-time dynamics of quarkonia, which becomes even 
more pregnant when the medium is rapidly expanding 
and cooling down, as is the case of the fireball in nucleus- 
nucleus collisions. 

For that purpose, it is interesting to consider alter- 
native modelings of the influence of the medium on the 
embedded QQ pair. Noting that the latter is a "small 
system" , then the surrounding medium can be seen as 
a "reservoir" , which can exchange energy and momen- 
tum with the small system without being noticeably af- 
fected. This is analogous to the paradigm setup for quan- 
tum dissipative systems [l2j], which suggests to view a 
quarkonium in a medium as such an open quantum sys- 
templEi. 

Accordingly, it becomes natural to study the dynami- 
cal evolution oi QQ pairs in a medium with the help of 
the techniques developed in the context of quantum dis- 



sipative systems. In a forthcoming paper [15[, we shall 
consider a wavefunction-based approach to obtain the dy- 
namics of quarkonia. Here, we use the master-equation 
formalism, and derive the time evolution of the popu- 
lations oi QQ states as well as of the momentum 
distribution of the quarkonia. For the latter, we show 
that it satisfies a Fokker-Planck equation, with trans- 
port coefficients fixed by the microscopic model. 



Various kinetic frameworks for the dynamics of quarko- 
nia in a medium have been considered in the litera- 
ture, based on the Boltzmann [l6l - fl8| , Fokker-Planck or 
Langevin equations [l9|, H(| or rate equations [2l|, [22j . 
to model the destruction and (re-) formation of bound 
QQ states in phenomenology-oriented studies. Here the 
open-quantum-system approach to heavy quarkonia in a 
medium provides a natural underlying microscopic de- 
scription that leads to such a kinetic model. 

In Sect.|lTl we introduce the general theoretical frame- 
work that we shall afterwards apply to obtain a micro- 
scopic description of the in-medium evolution of quarko- 
nia. Section IIHI introduces our model for the heavy 
quarkonia, the medium, and their interaction between 
them. Since we intend the present study to be of 
exploratory nature, we deliberately adopt a simplified 
model, instead of more realistic ones. Our results for 
the evolution of both the internal and external degrees of 
freedom, namely the populations of the various states 
and the QQ-vair momentum distribution function re- 
spectively, are presented in Sect. [IVJ Finally, we discuss 
our model together with the underlying assumptions and 
our results in Sect. |V1 where we also consider how these 
results might be modified in an evolving medium. 



II. EVOLUTION EQUATIONS 

In this section, we briefly review the master-equation 
description of quantum dissipative systems for the sake 
of self-containedness (a longer presentation can be found 
e.g. in Ref. After introducing in Sect. Ill Al the 

generic setup and its description, we present the equa- 
tions that govern the evolution of the dissipative system, 
starting with its internal degrees of freedom (Sect. Ill B[) 
and then turning to the external ones (Sect. QTCl) . 



A. Small quantum system coupled to a reservoir 

Generically, the setup for a quantum dissipative sys- 
tem consists of a (small) system S coupled to another 
quantum system 1Z, called environment — or reservoir, if 



2 



it has infinitely many degrees of freedom, as we shall as- 
sume from now on. The total system S + 1Z is assumed to 
be closed. It is then described by a Hamiltonian, taken 
to be of the form 

H = H S + H n + V (1) 

where H$ denotes the free Hamiltonian of the small sys- 
tem (in the absence of the environment), H-jz is the free 
Hamiltonian of the reservoir, and V describes the inter- 
action between system and environment. 

Hereafter, we shall model the reservoir as a set of har- 
monic oscillators, labeled by a subscript A, whose proper 
frequencies u\ span a large continuum, encompassing the 
Bohr frequencies of the free Hamiltonian Hs ■ Let p n de- 
note the density operator of the free reservoir. 

For the interaction term in the Hamiltonian, we con- 
sider a coupling of the form 

V = SR with R = J2(g x a x +g* x a{), (2) 

A 

where S acts on S only, while a\ and a\ are the anni- 
hilation and creation operators for oscillator A, and g\ 
measures the corresponding coupling. 

For a large reservoir, the autocorrelation function 
(R(t) R(t — r)) takes non-negligible values only in a small 
interval around r = 0, of typical size r c . 

Evolution equations for quantities pertaining to the 
small system are conveniently obtained by introducing 
first the density operator p of the total system, whose 
evolution is then governed by Heisenberg's equation with 
the Hamiltonian (jXJ) . Iterating the latter (in the Dirac 
interaction picture) and performing a partial trace over 
the reservoir degrees of freedom, one finds an exact, yet 
non-local in time, evolution equation for the "reduced" 
density operator 

/(t)=T> w (p(t)). (3) 

To obtain more tractable equations, a few simplifying 
hypotheses are needed. The first one consists of assuming 
that the total density operator is factorizable at any time: 

p{t) ~p s {t)®p n , Vt. (4) 

This amounts on the one hand to neglecting correlations 
between the small system and the reservoir beyond a 
certain order in the interaction term — typically, beyond 
second order. On the other hand, keeping the free den- 
sity operator p n in the presence of an interaction with 
S amounts to assuming that the latter does not mod- 
ify the reservoir properties, which is reasonable for the 
application we have in mind. 

The second assumption is that the typical time scale 
for the evolution of the small system should be much 
larger than the typical time scale r c of the reservoir fluc- 
tuations. 

Under these two hypotheses, one can derive a first- 
order differential "master" equation for the reduced den- 
sity operator p s {t) (23|. 



B. Evolution of the internal degrees of freedom 

In a first step, one can focus on the evolution of the 
internal degrees of freedom of the small system, momen- 
tarily leaving aside its motion. 

Let \i), \j) . . . denote the eigenstates of the free Hamil- 
tonian Hg, with Ei, Ej . . . the corresponding energies. 
In the basis spanned by \i){i\, the elements of 

the reduced density operator obey a set of coupled first- 
order differential equations with constant coefficients. 

For our purpose, it is sufficient to consider the diag- 
onal elements pf it corresponding to the populations of 
the energy eigenstates. These populations satisfy cou- 
pled Einstein equations of the form 

M( t) = -^Ti^kpUt) + J2 T ^iPkM (5) 

with transition rates given (when considering the master 
equation up to second order in the interaction term) by 
Fermi's golden rule. For Ek > Ei and ujki = (£* — Ei)/% 
the corresponding Bohr frequency, one easily finds 

r fe ^ = |^^((n A ) + l)|(i;l A |y|fc;0)| 2 5(w A -a; fei ), (6a) 
n x 

r '^fe = S X>A> I <^ !a) -w«), (6b) 

where (n\) denotes the average number of excitations in 
mode A. These rates obviously correspond to emission 
(with the +1 term accounting for spontaneous emission) 
and absorption, respectively. 

To exemplify the behavior of populations described by 
Eqs. (IS])-©, we show in Fig.[T]the time dependence of the 
populations of a 4-level system S, initially in its ground 
state, in contact with a thermal bath at temperature T. 
In that case, (n\) is given by the average occupation 




time (arbitrary unit) 



FIG. 1. Time dependence of the populations of the states of 
a 4-level system coupled to a thermal bath. The straight lines 
correspond to the equilibrium values at the bath temperature. 
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number for Bose-Einstein statistics. After a transient 
regime, the populations reach stationary values, which 
are simply proportional to the corresponding Boltzmann 
factors: 



rii / cq. 



Eh — Ei 
k R T 



We thus find that the internal degrees of freedom of the 
small system equilibrate at temperature T. 



C. Evolution of the external degrees of freedom 

If we also consider the motion of the center of mass 
of the small system, then momentum transfers due to 
emission or absorption of excitations from the reservoir 
now play a role. In addition, the dependence of the in- 
teraction term ([2]) on the position X of the small system 
should now be taken into account. For the case we shall 
be interested in later on, this can be done by replacing 
a x by a A e lk * x with |k A | = u\/c. 

The main difference with the static case is that the 
eigenstates of the free Hamiltonian H$ are states with 
not only different internal quantum numbers, but also 
different momenta, corresponding to the P 2 /2M$ part 
of Hs, with Ms the mass of the small system and P 
its total-momentum operator. Consequently, the eigen- 
states should be labeled with the eigenvalue p besides 
the "internal" label i. 

Introducing now the shorthand notation 7T; iP = pf- t pp 
for the diagonal terms of the reduced density matrix — 
which can be viewed as momentum distributions when 
they are considered as functions of p — the evolution 
equations flSJ) become 



d7T 



At 



£(*) 



■El 



rfc,q-M,p7rk,q(*) — Ti -p ^k,q^i,p (t) , (7) 



with rates again given by Fermi's golden rule. In evalu- 
ating the latter, it is convenient to consider the position- 
dependent part of the interaction Hamiltonian apart 
from the rest. This part, convoluted with the position- 
dependent part of the Hs eigenstates — namely plane 
waves — gives rise to a momentum-conservation enforcing 
term <5 Piq _/jk A in the matrix element (k, q; 1 p; 1 A ). 
Besides, one should also include the kinetic energy con- 
tributions. All in all, one finds for the transition rates 
between two levels with Ek > Ei 

r*, P ^, q = ^5>i A }|<fc;0|V|i;l A }| 2 



A 



x d Pi q-Hh x 5(u!\ - Lu kl - t D + £ R ), (8a) 



A 

x <5 p , q _nk A <S(wA - u ki - £d - £r), (8b) 



where V denotes the position-independent part of the 
interaction term. £d and £r are the frequency shifts due 
to the Doppler effect and the recoil effect, respectively: 



c A • p UJ\ 



M s 



2 Afjc 2 ' 



2M S ' 



(8c) 



where the Doppler effect has been considered up to sec- 
ond order. Inserting the rates (JHJ) into Eqs. (0, the sums 
over q disappear thanks to the momentum-conservation 
condition. 

Further analytical progress with the population evolu- 
tion equations requires additional assumptions, namely 
first that the frequency shifts £d, £r remain much smaller 
than the typical width Aw of the bath spectral distribu- 
tion; and, secondly, that the momentum transfer hk\ be 
much smaller than the width Ap of the populations 7Ti, P , 
viewed as momentum distributions. 

Let the sum of the populations 7Ti. p over all inter- 
nal states % be denoted by 7r(p), which then represents 
the momentum distribution function of the small sys- 
tem, irrespective of its internal state. As shown in Ap- 
pendix IA 31 under the assumptions mentioned above the 
rate of evolution for 7r(p, f) is much slower than the in- 
dividual rates d^p/di, and one can show that 7r(p,t) is 
governed by 



at 



VdV p ■ [p7r(p, t)] + k A p7 r(p, t), (9) 



i.e. an equation of the Fokker-Planck type, in momentum 
space, with constant coefficients t\d and k. The former 
describes the damping rate of both the average momen- 
tum and (up to a factor of 2) the variance of the mo- 
mentum distribution, while the latter characterizes the 
growth of this variance. When the reservoir in which 
the small system evolves is a thermal bath, both coeffi- 
cients are related to each other through the fluctuation- 
dissipation relation 



Ms 



T] D k B T, 



(10) 



which shows that in the stationary regime, the small sys- 
tem has thermalized at the bath temperature T. 



III. MODELING QUARKONIUM AS A 
QUANTUM DISSIPATIVE SYSTEM 

As stated in the introduction, our goal in the present 
paper is not to propose a full treatment of the dynamics 
of heavy quarkonia in a thermalized medium based on the 
most refined existing models for quarkonia and their in- 
teraction with the medium. Our purpose is rather to ex- 
plore possible new qualitative phenomena, which emerge 
when the point of view on the problem is shifted from 
the usual approaches to the description as a quantum 
dissipative system [T3 |. 

For that reason, the models we introduce hereafter 
for the quarkonia (Sect. MI Bp and their coupling to the 
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medium (Sect. IIII C|) will be quite simplified, yet not 
unrealistically. For the "reservoir" with which the QQ 
states interact, we consider two possibilities: either a 
thermal bath, or a peaked distribution (Sect. IIII A|) . 



A. Medium as a reservoir 

In nuclear collisions at sufficiently high energies, the 
medium which is created and which is probed by heavy 
quarkonia should be deconfined, and thus consist of 
quarks and gluons as relevant degrees of freedom. For the 
sake of simplicity, we consider a medium made of pure 
glue, and forget the constantly created quark-antiquark 
pairs, which would not affect the qualitative features of 
our description. We assume that this gluon plasma is 
unpolarized, isotropic and homogeneous. 

This plasma can then be modeled — for example by 
quantizing the SU (3) gauge fields canonically in the 
Weyl gauge — as a set of oscillators. Within the master- 
equation formalism, the only characteristic we need is the 
mean number of excitations (n\) for each mode A, see the 
transition rates © and ©• 

In the following, we shall make use of two different 
kinds of gluon bath. First, we shall consider a thermal 
bath, i.e. a reservoir in thermal equilibrium, at a temper- 
ature T. The associated density operator reads 



c -HK/k B T 

Tr(e- ff «/^ T ) 



(11) 



in the absence of the small system. The corresponding 
(n\) is given by the usual Bose-Einstein distribution. 

The second model of reservoir we shall employ consists 
of assuming a Gaussian distribution 



(n\) oc exp 



hu>\ — E 
2AE 



(12) 



peaked around some variable value E, for the average 
number of excitations. Although this is of less immedi- 
ate relevance for the phenomenology of ultrarelativistic 
heavy-ion collisions than the thermal bath, yet it will al- 
low us to illustrate some features of our description. This 
second model will be referred to as the "Gaussian bath" . 



B. QQ states 

Strictly speaking, to implement the master-equation 
formalism described in Sect. [TT] one only needs matrix 
elements for reservoir-induced transitions between states 
of the small system. Specifying the states themselves and 
the transition-inducing interaction is not necessary. 

Accordingly, for the study of the dynamics of heavy 
quarkonia in a gluon plasma, one should identify all 
single- or multi-gluon processes that change the state 
of a QQ pair — be it a transition between two different 



bound states [24j , gluon- induced dissociation [25j, |26| , or 
the possible recombination of a quark and an antiquark 
into a bound state [27j — , and consider the corresponding 
matrix elements. Such an exhaustive procedure is cer- 
tainly desirable for making quantitative predictions that 
can meaningfully be compared to experimental results. 
Here we remain at an exploratory level, and search the 
qualitative behaviors of quarkonia in a medium. To make 
amend for our not using the most accurate set of matrix 
elements, we do not restrict ourselves to postulating such 
a set, but we shall start from scratch, i.e. from a model 
of quarkonia in the vacuum, and of their interaction with 
the gluon plasma introduced above. 

For the purpose of identifying new phenomena, the 
bottomonium system, with its larger number of bound 
states likely to survive above the deconfinement temper- 
ature seems to us more promising than the charmo- 
nia. A further advantage of bottomonia, is that (in vac- 
uum) they can reasonably be described as bound energy 
eigenstates, with simple wavefunctions, of a one-gluon- 
exchangc Coulomb potential 



V Q o(r) 



-C F 



a s hc 



(13) 



with Cf =4/3 the usual color factor and a s the dimen- 
sionless (running) coupling constant, here a s ~ 0.25. For 
charmonia, this would be a less satisfactory description. 

There are several drawbacks to our modeling bottomo- 
nia as bb pairs bound by a Coulombic potential. First, the 
eigenstates of potential (fTS"]) come in degenerate subsets, 
while this degeneracy — which prevents direct transitions 
between degenerate states — is lifted in the correspond- 
ing bottomonia. To allow such direct medium-induced 
transitions, we lift the degeneracy by hand, and give the 
states their vacuum masses (28j. 

Another issue is that not every known bottomonium 
has been assigned all its quantum numbers, and some 
states (e.g. in the ID-quintuplet) have not yet been 
identified experimentally. To deal with these "missing 
states" , we retain the degeneracy of states within P-wave 
triplets and D-wave quintuplets, even when they are dif- 
ferentiated experimentally. 

Thirdly, while potential (|13D admits an infinity of 
bound states, only a handful of bottomonia are actu- 
ally stable against the strong interaction. And last, even 
though the scattering states of the Coulomb potential 
are known, yet we found it disturbing to use them to 
describe transition processes (dissociation or recombina- 
tion) between bound bottomonia and free (anti)quarks, 
given that the latter do not exist in the vacuum. To cope 
with both these problems, we made a drastic assump- 
tion, namely that the bound eigenstates of potential (Ti~3"|) 
above and inclusive the 45-level stand for unbound bb 
states. Additionally, we forbid by hand transitions from 
such "unbound" states back to bound ones. 

The resulting spectroscopy of states we consider, with 
the transitions between them which we shall detail in 
Sect. IIII Cl are displayed in Fig. [21 Note that this slightly 
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FIG. 2. Scheme of the spectroscopy of 66 states and their 
transitions used in the calculations. 



differs from the spectroscopy we used in [14j , inasmuch as 
we have now added the D-wave states, which will impact 
our results due to their large overlap with the P-wave 
states. To estimate the error on our results, we also add 
one further level of (unbound) states, to which the bound 
levels can transition, but which cannot transition back. 



of gluon A, while d (resp. r) denotes the dipole (resp. 
radius) operator for the QQ pair. 

Such an interaction term induces, to first order, transi- 
tions between QQ states with different orbital quantum 
numbers only, as represented in Fig. [2] 

Note that the dipolar coupling ([XI)) actually rests on 
the assumption that both quark and antiquark in the pair 
see the same chromoelectric field. That is, we implic- 
itly assume a large wavelength in computing the rates. 
While this holds for the gluons that induce transitions 
between bound states, yet it might not be granted for glu- 
ons that would dissociate the ground quarkonium state — 
for which one could instead use the rate computed in 
Ref. fUl, which we did not do. 



IV. RESULTS: EVOLUTION OF HEAVY 
QUARKONIA IN A THERMAL MEDIUM 

Within the framework of the model we have introduced 
in the previous section, we can now turn to solving the 
evolution equations for the populations of bottomonia 
in a thermal medium at temperature T. Following the 
same order as in Sect. [Til we first present results for the 
evolution of the internal degrees of freedom (Sect. ITVAl) . 
then for the dynamics of the quarkonium center of mass 
(Sect. HVBp . In Appendix [Bj we present results for the 
evolution in a Gaussian bath. 



C. Quarkonium-plasma interaction 

Eventually, we need to specify the interaction between 
a QQ pair and the gluon plasma. In this work we restrict 
ourselves to considering dipolar coupling, which induces 
vector transitions in the QQ system^ This amounts to 
considering the coupling of the QQ pair to the dipolar 
part of the chromoelectric field of the gluons, which for 
an unpolarized plasma yields the interaction term 

V = — d E = -\y/C F a. s hcr ■ ^2 y ^ffi £ a(qa - a{), 

(14) 

with L the size of the box in which the chromoelectric 
field is quantized (which also appears in the normaliza- 
tion of the bb eigenstates) and e\ the polarization vector 



One might worry that singlc-gluon interactions induce transitions 
from color singlet to color octet QQ states. While this is certainly 
true, yet it should be kept in mind that the model of quarkonia 
as pure bound QQ states is only approximate: taking account 
the sea, an improved picture for a quarkonium is rather 

l(QQ)> = ^ Q q\QQ) + 4> Q Q 3 \QQg) + ^ Q Q qq \QQqq) + ■■■ . 

where the QQ pair in \QQg), \QQqq)... can be in the octet repre- 
sentation, i.e. each quarkonium actually contains some admixture 
of color octet QQ. 



A. Internal evolution of a static quarkonium 

Inspecting the set of equations (O that govern the be- 
havior of the populations when medium-induced emission 
and absorption processes are taken into account, a first 
result strikes the eye, even before solving the equations. 

Let the populations p®® be combined into a vector p and 
the system ([5]) be rewritten as 

with Utz the time-evolution operator for the populations. 
In the vacuum, the matrix representation of U-n m the 
basis of the energy eigenstates of the QQ system (ordered 
in increasing energies) is triangular. This is no longer 
the case in the presence of a medium. Consequently, the 
eigenvalues and eigenstates of Uu are not the same in a 
medium as in its absence^ Physically this implies that 
the higher-energy QQ states do not evolve independently 
from the more bound ones, as in the vacuum: because of 
medium-induced transitions, the former become coupled 
to the latter. As a consequence, past a transient regime, 
the populations of all states evolve with the same time 
scale. 



2 For a proof that Uji is indeed diagonalizable, see Ref. 
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FIG. 3. Evolution of bottomonium populations in a thermal 
bath at T = 5T C , with T c = 170 MeV. 

This holds irrespective of the assumed matrix elements 
between the QQ states — provided every state is coupled 
to each other, at least indirectly. However, the value of 
the time scale is model-dependent and depends on the 
matrix elements and on the bath properties. In Fig. [3] 
we show the time evolution of bb states, modeled as in 
Sect. HID in a thermal bath at T = 5T C , where the as- 
sumed initial condition consists of having all pairs in 
the ground state T(15) at t = 0. The curves do not 
change significantly if we include one further level of un- 
bound states (not shown). After the first fm/c or so, 
one reaches a quasi-equilibrated regime where the popu- 
lations of all vacuum bound states decay with a charac- 
teristic time scale of 1.5 fm/c, while their ratios remain 
stationary!! This result should be contrasted with the 
sequential-melting picture, where states would either be 
there, or totally melted according to the medium tem- 
perature, but cannot be regenerated through transitions 
from other states. 

In a thermal bath at 2T C , we find 8 fm/c for the time 
scale of bottomonium evolution: as could be anticipated, 
the time scale decreases with rising temperature. 

Focusing on the quasi-equilibrium population ratios, 
we show in Fig. [4] their dependence on the temperature 
of the plasma. These ratios differ significantly from their 
values for thermally equilibrated bottomonia, as would 
be expected in the framework of a statistical model [3(| ■ 
This difference can easily be traced back to our forbidding 
transitions from unbound states to bound ones, so that 
the detailed balance condition 

e~ E ^ T = T k ^ e- E ^ k - T Vi, k, (15) 

which guarantees the existence of an equilibrium with 
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FIG. 4. Temperature dependence of the ratios of bottomo- 
nium populations. Symbols: quasi-equilibrium ratios within 
the master-equation formalism; full lines: ratios in a ther- 
mally equilibrated system. 

populations proportional to the respective Boltzmann 
factors (see Sect. Ill 51) . does not hold here. 

One could argue that we have put this deviation from 
thermal equilibrium at long times by hand, by forbidding 
some of the emission transitions. This is true, but ulti- 
mately due to our over-simplified modeling of unbound 
states. As long as only a finite number of them is explic- 
itly included, together with the back transitions, then 
equilibrium is reached after some finite time, which in- 
creases very rapidly with the number of states^ When 
unbound states form a continuum, equilibrium is reached 
infinitely late, which is what we have modeled by setting 
some transition rates to 0. 



B. Evolution of the external degrees of freedom 

Now that we have understood the internal dynamical 
evolution of the quarkonia, we can turn to investigating 
the evolution of the external degrees of freedom, and es- 
pecially of the momentum distribution. 

More precisely, we wish to consider the dynamics of the 
"momentum distribution of bound quarkonia" 7r(p, t), 
denned as the sum over all bound levels of the popu- 
lations 7Tj lP (t) = /0^pp(i) ; where the density matrix is 
taken in the basis of the (vacuum) energy eigenstates. 
This distribution can evolve under the influence of two 
different effects induced by the plasma. 

First, the QQ bound state can be dissociated, i.e. it is 
"lost" from the populations that enter 7r(p, t), which thus 
decreases with time. More precisely (see Appendix IA 2j) , 



3 The time scale reported in our previous work \li was larger 
because there we had not considered the D-wave states. 



4 The attained equilibrium might differ from the thermal values 
when condition H15B is not fulfilled. 
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When transitions between the various quarkonium 
states are allowed, then in the presence of a medium 
at fixed temperature, after a transient regime a sta- 
tionary stage is reached, in which the populations 
of all states evolve together, as illustrated by Fig. [3] 

The ratios of these quasi-equilibrated populations 
in the stationary regime differ from the ratios in a 
statistical model for quarkonia in thermal equilib- 
rium with the plasma, see Fig. 01 
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FIG. 5. Drift (full line) and diffusion (dashed line) coefficients 
of the Fokker-Planck equation ([9]) for the bb system immersed 
in a thermal bath at temperature T. 

7r(p, t) decays exponentially, with a rate that depends on 
p. The latter point is easily understandable, inasmuch 
as the energy of the gluons that interact with the moving 
QQ pair is Doppler-shifted, so that pairs with different 
momentum do not see the same gluon spectrum. 

On the other hand, the gluons can also induce inter- 
nal transitions between bound states of the QQ pair. In 
that case, the master-equation formalism predicts that, 
at least at small momentum |p| and for small momen- 
tum transfers 7i|k|, the rate of change of 7r(p,i) is sig- 
nificantly slower than the rates of the individual pop- 
ulations. Furthermore, in this regime the momentum 
distribution obeys the Fokker-Planck equation ©. We 
display in Fig. [S] the temperature dependence of the drift 
and diffusion coefficients in this equation, computed for 
the same bottomonium system as in Sect. lFVAl in a ther- 
mal bath0 As can be checked, these coefficients satisfy 
the fluctuation-dissipation relation (ITUl) . While this hints 
at the equilibration of the external degrees of freedom of 
the bottomonia in a thermal plasma, yet one should keep 
in mind that the relevant time scale ~ r]^ 1 might actually 
be significantly larger than the time scale for bottomo- 
nium dissociation. 



V. DISCUSSION 

In this paper, we have applied the master-equation for- 
malism to describe the evolution of heavy quarkonia in 
a gluon plasma, in complete analogy with the evolution 
of a small quantum system in contact with a reservoir. 
Independently of the model we used for the quarkonia, 
the plasma, and their interaction, several generic features 
emerge: 



5 For such a bath, the condition on momentum underlying the 
derivation of the Fokker-Planck equation amounts to requiring 
that the bottomonia be non-relativistic. 



The momentum distribution of bound quarkonia, 
considered irrespective of the internal state, satis- 
fies a Fokker-Planck equation, at least in the non- 
relativistic regime. 

A further expected behavior — which we have not inves- 
tigated in the present work, but manifests itself when 
studying the non-diagonal elements of the density ma- 
trix [23J — is that the interaction with the medium shifts 
the energy levels of the small system with respect to their 
vacuum values. 

For the sake of illustration, we considered a simpli- 
fied model for the heavy quarkonia, and more partic- 
ularly for bottomonia rather than charmonia, and for 
the medium-induced transitions. Despite the rudimen- 
tary character of these models, the numerical values that 
come out for the characteristic time scale of the evolution 
of bottomonium populations and for the drift coefficient 
in the Fokker-Planck equation, including their depen- 
dence on temperature, are actually very similar to the 
values derived in more elaborate models for the T(IS') 
lifetime [3l| and for transport coefficients in the hard- 
thermal- loop approach (32j | . This is an encouraging find- 
ing, that shows the potential of the approach. 

Here we wish to emphasize again that the starting 
point for the implementation of the formalism is either 
transition rates or, if one wants to start from scratch, a 
description of the quarkonia in vacuum and of the inter- 
action with the medium. It is certainly tempting to use 
an in- medium quark-antiquark potential [33j-|38[ . Yet one 
should not forget that such a potential already accounts 
for part of the coupling to the plasma, which then has to 
be subtracted out in a consistent manner from the other 
ingredients of the model, to avoid double counting. For 
instance, one can qualitatively expect that the screening 
of the potential in an in-medium potential picture corre- 
sponds to an increase of the transition rates from bound 
to unbound states in the formalism of the present paper. 

Reformulating this differently, the master-equation for- 
malism provides an evolution equation for the reduced 
density operator describing the quarkonia. Under the as- 
sumptions that make it local in time, this equation might 
be equivalent to a Heisenberg equation for pQQ involving 
a Hamiltonian with an effective potential, which incorpo- 
rates the influence of the plasma, irrespective of whether 
the latter is in thermal equilibrium or not. Since we used 
a simplified vacuum quark-antiquark potential, we have 
not attempted to extract some corresponding effective 
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t (fm/c) 

FIG. 6. Evolution of bottomonium populations in a thermal 
bath with evolving temperature. 

potentially Even then, it is clear that this in-medium 
potential would include an imaginary part, to account 
for the non-unitarity of the evolution of p®Q ■ 

For a future application to a more precise description 
of heavy quarkonia and their interaction with the fireball 
created in ultrarelativistic nucleus-nucleus collisions, one 
should discuss two aspects, namely the validity of the as- 
sumptions underlying the master-equation formalism and 
the relevance of the features listed above in the context 
of interest. 

The main hypotheses behind the derivation of the mas- 
ter equation are twofold. First, the formalism holds pro- 
vided the characteristic time scale of the medium fluctu- 
ations is much smaller than the time scale of the small 
system dynamics. This ensures that the possible corre- 
lations between medium and small system are continu- 
ously washed out, so that the evolution of the latter is 
Markovian. Given the size of the medium — equilibrating 
"parton" gas, quark-gluon plasma or hot hadron gas — 
created in high-energy heavy-ion collisions, this point is 
warranted. The second hypothesis, namely that of a 
"weak" coupling, which underlies the use of transition 
rates given by Fermi's golden rule, is actually less cru- 
cial. In a forthcoming study [l5j], we shall introduce an 
alternative approach releasing this assumption; however, 
the coupling strength does not affect the qualitative re- 
sults summarized above. 

Even though the master-equation approach seems to 
be applicable, at least as a good approximation, yet for 
the evolving medium created in high-energy heavy-ion 
collisions, some of those results have to be reexamined. 
Thus, the rate of evolution of the medium might be com- 
parable to the rates of quarkonium-plasma interaction, 
and prevent the equilibration of the internal and exter- 



This would amount to performing a Kraus decomposition [39l of 
the mapping from pQQ(0) to p ( 3'3(t). 



nal degrees of freedom of the quarkonia. As an example, 
we present in Fig. [B] the evolution of the populations of 
bottomonia in a deconfined plasma whose temperature 
decreases with time as found at the center of the interac- 
tion region in hydrodynamical simulations of Pb-Pb col- 
lisions at LHC energies [4(j . One can see that the various 
bound states do not evolve similarly, so that the ratios of 
their populations do not remain constant, in contrast to 
point 1 above. Similarly, the description of the evolution 
of the momentum distribution of bound states through a 
Fokker-Planck equation, which relies on the equilibration 
of the internal degrees of freedom, does not hold either if 
the fireball cools down too rapidly. 

This shows that results derived within a stationary pic- 
ture for the quarkonia and/or the medium might actually 
not hold when real-time evolution is taken into account. 
The master-equation formalism, and other descriptions 
of dissipative quantum systems, can accommodate such 
a real-time evolution, since the assumptions made about 
the medium are rather minimal. Here we have shown that 
such an approach is possible, yet an accurate modeling 
of the dynamics of heavy quarkonia in ultrarelativistic 
nuclear collisions deserves further investigation. 



ACKNOWLEDGMENTS 

We thank Nirupam Dutta and Helmut Satz for helpful 
discussions and suggestions. C. G. acknowledges support 
form the Deutsche Forschungsgemeinschaft under grant 
GRK 881. 



Appendix A: Motion equation for quarkonia in a 
medium 

In this appendix, we detail the derivation of the 
Fokker-Planck equation that describes the evolution in a 
medium of the momentum distribution 7r(p,t) of bound 
quarkonium states. 

As stated in Sect. IIIH within our model one has to dis- 
tinguish between transitions between bound QQ states, 
which can take place in both directions, and transitions 
between a bound state and an unbound one, which can 
be dissociation processes only. 

For the sake of brevity, we hereafter consider two 
different systems — coupled to reservoirs — with each two 
non-degenerate levels E a < E^. Let 7r a!P , 7rf, iP denote 
the diagonal elements of their respective reduced den- 
sity matrices in the energy eigenstate basis. In the first 
system ("system I"), both levels correspond to "bound 
states" that can transition to each other. We then call 
"momentum distribution of the bound states" the sum 

7r r (p, t) = 7T Q , p (t) + Wb,p(t). 

In opposition, for system II, only excitations from the 
lower to the higher level are allowed, while transitions 
back are forbidden, so that this constitutes an open sys- 
tem. The "momentum distribution of the bound states" 
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is then ir ll (p,t) = w a>p {t). Generalizing the calculation ward and amounts to combining the two behaviors that 
to more complicated spectroscopies is then straightfor- we encounter below. 



1. Evolution equations 

For those systems, the evolution equations ([7]) with the transition rates ([5]) read 



dir. 



2tt 



dt 



(Ala) 



(Alb) 



where LOb a is the Bohr frequency of the transition, £d and £r are given by equation (|8cl) . while rj = 1 for system I, 
7] = for system II Q 

Introducing the quantity 



fab = 77 5Z ((MI^I^Ia)! 2 ^^ -w ba ) 



and a continuum representation, the evolution equations (lAip can be rewritten as 



d"7T„ 



dt 

d7Tfe ; p t 

At 



(t)=f a b I I 7— S{uj - u ba — £ D - £ R ) (rj [(n(w)) + lj^p+skC*) - («(<*>)} 7r 0>p (i)) ; 



dft 



(A2) 

(A3a) 
(A3b) 



I 

where the modulus |k| = u/c is fixed by the respective with 
Dirac distribution — noting that £d and £r actually also 
depend on u) — , while dJ7 is the elementary solid angle 
around the direction in which k is pointing. 

These equations clearly show that the motion of the 

(A5) 

small system center of mass shifts the frequency of the where the dependence on p in the right-hand side is hid- 
absorbed or emitted bath excitations away from the Bohr den in theDoppler frequency shift £ D . The solution to 



a— >b,p 



-r{n(u)} [ 7^ 5(w-u ba -£ D -£ R ), 



frequency. These processes also lead to an increase or a 
decrease of the momentum of the small system by Kk. 



2. Open system: exact solution of the evolution 
equation 

When rj = 0, that is for system II, equations (IA3[) be- 
come simpler, and in particular the evolution of the pop- 



equation (|A4[) is trivial. 



3. Closed system: perturbative expansion of the 
evolution equation 

For 77 ^ 0, in particular for system I, further analytical 
progress with equations (|A3I) necessitates extra condi- 
tions on the size of the frequency shifts £d, £r and of 



ulation of the "bound state" decouples from that of the the momentum transferred in an absorption or emission 
higher level. Equation (|A3ap , which describes the evolu- 
tion of the bound-state momentum distribution tt 11 = n a , 
now reads 

dTr 11 

— (p,i) = -r^pTr 11 ^,*), (A4) 



process, namely £d, 6? "C Aw, the width of the reservoir 
spectrum, and fik <C Ap, the width of the momentum 
distribution. To ensure that these assumptions hold, it 
is sufficient that the two parameters 

ei = — and e 2 = . A6) 

More generally, a factor 77 ^ 1 might account for non-equilibrated be much Smaller than unity, since this implies automati- 



up and down transition rates. 



cally £r/Alj - Sie 2 < 1. 
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Under these assumptions, we can Taylor-expand up to 
second order the Dirac distribution 

5(oj ba - u+£ D ±£ R )~8(u ba - w) + (£d±£r) S'(u ba - u) 



5"(u) ba -cj) (A7) 



and the momentum distributions 7r a p , 7Tb !P 

h 2 d 2 ir 

^^kW - 7T P (t) ± Hk • V7r p (i) + — J2 fcfcj«r»r-(*) 



2 ""' dpidp, 

1,3 



(A8) 

in equations (|A3I) , and deduce simplified evolution equa- 
tions for the bound-state momentum distribution ir l by 
identifying the factors of the various powers of e\ and £2. 
The zeroth-order terms, which amount to neglecting 



the momentum transfer and the frequency shifts, are triv- 
ial and express the global conservation of the population 
of the system: d7r I (p)/dt = 0. In turn, the linear terms 
in £1, £2 are automatically proportional to k, and thus 
yield a vanishing contribution when averaged over all di- 
rections for k. 



At quadratic order, there are terms in e\ — from the 
second-order term in expansion (|A8I) — and in £i£2 — from 
the term linear in £r in (IA7|) and from the product of the 
first-order terms of both expansions. The terms in e| can- 
cel out when summing equations (|A3ap and (|A3b]). All 
in all, after performing the integrations by part necessary 
to get rid of the 6' terms, followed by the straightforward 
integrations over u and the direction of k, one obtains 



dTT 1 , s ft 



ba 



f[(n(uj ba )) + l] Air b . p (t) + (n(u) ba )) A7r a , p (i) 



fa>L d ( n ), 



- 5hu} ba 

1 abTTTT -V 



3M 5 C 2 



p([{n(ui ba )) + l]ir bjP (t) - (n(uj ba )) 7r 0)P (f) 



(A9) 



The problem with this evolution equation is that it 
still contains the individual internal populations of the 
various internal states of the small system, and not only 
the total populations. An equation involving only the 
latter can be derived, provided the internal populations 
7r a,p 7 Tb, p remain in fixed ratios, i.e. when the internal 
degrees of freedom are equilibrated. 

Now, the evolution rates d7r iiP /di themselves are of 
order in s\, £2, much larger than the evolution rate 
for 7r(p). That is, one may assume that the internal 
degrees of freedom reach stationary values on a much 
smaller time scale than the typical scale for the evolu- 
tion of the momentum of the system. Inspecting equa- 
tions (|A3I) with vanishing left-hand sides and considering 
only the leading terms in the Taylor expansions (|A7[) and 
one checks that the prescriptions 

1 + (n{u ba )) 



7T5, p (t) 



1 + 2(n(uba)) 

(n(u ba )) 
l + 2{n(u ba )) 



7r(p,i) 



are stationary solutions the evolution equations at lead- 
ing order. Inserting them in equation (|A9|) . one obtains 



dTT 1 



dt 



(p,i) = /tA7T I (p,t)+7 ?D V- [pTr^p,*)]. (AlOa) 



(n{u ba ))}(n(ui ba )) H u>l a ~ 



1 + 2{n{u) ba )) 



3c 2 



r a b 



(AlOb) 



and 



TvJl a d(n) 



1 



where we have set 



.>, , > , (Wba)fab- (AlOc) 

2(n(w 6a )) 3M s c 2 du 

Equation (IA10a[) is an equation of the Fokker-Planck 
type, with constant diffusion coefficient k and drift coef- 
ficient r\Y) . Note that if one pushes the Taylor expansions 
to the next order, then the extra terms can be reexpressed 
as a momentum dependence of k and r\u- 

In the case when the reservoir with which the small 
system is in contact is a thermal bath at temperature T, 
then (n(w)) is given by the Bose-Einstein distribution 
function, and one easily checks that the diffusion and 
drift coefficients (|A10bj) - (|A10cl) satisfy k = M s k B Tr] D . 

Let us end this appendix with some remarks on the as- 
sumptions underlying the derivation of equation ([AlOajl . 
The Taylor expansion in powers of £1 and £2 is also the 
ingredient behind Landau's derivation of the Fokker- 
Planck equation from the Boltzmann equation without 
mean field term [4l|. The smallness of £1 then amounts 
to considering that soft momentum exchanges with the 
bath play the major role in the evolution of the motion 
of the small system. Eventually, since the frequency shift 
due to the Dopplcr effect increases with the system mo- 
mentum, the requirement of a small £2 implies that the 
description might not hold at high momentum. 
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Appendix B: Evolution of heavy quarkonia in a 
Gaussian bath 

In this appendix, we wish to take advantage of the 
fact that the master-equation formalism can accommo- 
date various models of reservoirs, not only thermal baths, 
and demonstrate that our model of quarkonia in a gluon 
plasma naturally incorporates a behavior which is ex- 
pected on physical grounds. 

For that, we immerse bottomonia in the ground T(1S) 
state at t = in the Gaussian bath (fT2")l . keeping the 
same models for bb pairs and their interaction with the 
plasma as in Sects. IIIIHVI As average energy E of the 



bath excitations, we consider first E = 5T C , then 10T C , 
with a width AE = T c in both cases. The resulting bot- 
tomonium populations, as a function of time, are shown 
in Fig. El 

One finds two very different behaviors. For E = 5T C , 
the populations evolve similarly to the case of a thermal 
bath at T — 5T C , see Fig. [5] On the other hand, for 
E = 10T C there is almost no evolution over the same time 
interval. The physical interpretation of the latter finding 
is simple, namely that the gluons in the second bath are 
too energetic to view the bottomonia as a whole, and 
therefore cannot excite or dissociate efficiently, as was 
already found (for J ftp gluodissociation) in Ref. [26}. 
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